Unveiling underestimated species diversity within the Central American Coralsnake, a medically important complex of venomous taxa

Coralsnakes of the genus Micrurus are a diverse group of venomous snakes ranging from the southern United States to southern South America. Much uncertainty remains over the genus diversity, and understanding Micrurus systematics is of medical importance. In particular, the widespread Micrurus nigrocinctus spans from Mexico throughout Central America and into Colombia, with a number of described subspecies. This study provides new insights into the phylogenetic relationships within M. nigrocinctus by examining sequence data from a broad sampling of specimens from Mexico, Guatemala, Honduras, Nicaragua, Costa Rica, and Panama. The recovered phylogenetic relationships suggest that M. nigrocinctus is a species complex originating in the Pliocene and composed of at least three distinct species-level lineages. In addition, recovery of highly divergent clades supports the elevation of some currently recognized subspecies to the full species rank while others may require synonymization.

www.nature.com/scientificreports/ and/or directly in the field. When specimens were obtained in the field (see Supplementary Table S1) and/or when euthanasia was necessary, the protocol and procedures employed (specimen capture, handling, and preservation) were ethically reviewed and approved by The University of Texas at Arlington Institutional Animal Care and Use Committee (IACUC office) (Title of protocol: Phylogeny of New and Old World Venomous Snakes; Protocol number: A07.031 to ENS; or "Studies of Reptiles and Amphibians" A08.025). Euthanasia involved an intraperitoneal injection of Pentobarbital, 0.1 ml/kg, or 10% chloretone cardiac injection, depending on specific country substance regulations. Field samples were preserved by snap freezing, kept dry (sheds or road kills), or preserved in lysis buffer, 100% EtOH or RNAlater (Thermo Fisher Scientific). All methods were performed in accordance with the relevant guidelines and regulations. The reported study is in accordance with ARRIVE guidelines. In total, we include 34 newly sequenced specimens (Table S1). DNA extraction protocol follows 32 .
To infer the phylogenetic relationships and reconstruct the timing of evolutionary divergence, we sequenced the mitochondrial protein-coding genes NADH subunit 4 (ND4, 666 bp) and cytochrome b (Cytb, 711 bp). Templates were sequenced on both strands, and the complementary reads were used to resolve rare, ambiguous base calls in Sequencher v. 4.9. Sequences were aligned in Seaview v. 4.2.11 33 under ClustalW2 34 default settings. Nucleotide translation into proteins had no stop codons. Primers are listed in Supplementary Table S2. The GenBank sequences generated for this study are given in Supplementary Table S1.
Phylogenetic analyses. We used PartitionFinder v.2 35 to choose the optimal partitioning strategy under a greedy search scheme 36 . The most appropriate substitution models for the Bayesian Inference (BI) and RAxML 37 analyses were determined by the Bayesian information criterion (BIC). Phylogenetic inference was based on the ND4 and Cytb gene fragments (complete alignment length; 1377 bp). The best partition scheme for the ND4 and Cytb alignment supported a codon partition for all analyses (BI, RAxML) (see supplementary Table S3). We only used previously published sequences with known vouchers, locality, and confirmed identification.
Bayesian inference trees were built in MrBayes v. 3.2 38 under substitution models with codon-based partitions, default priors, and Markov chain settings, and random starting trees. Each run consisted of four chains of 20,000,000 generations, sampled every 20,000 generations, discarding 25% of the trees as burn-in. We used the software RAxML v7.0.4 37 to perform a maximum likelihood (ML) inference of the evolutionary relationships between our samples with default settings and codon partition (Supplementary Table S3). Tree inferences were performed through the CIPRES platform 39 . In both analyses, we assigned Micruroides euryxanthus as the outgroup 40 . Genetic p-distances of combined Cytb and ND4 were calculated in MEGA X 41 under a pairwise partial deletion option.
Fossilized birth death divergence times estimation. We estimated divergence times using the Sampled Ancestors (SA v2.0) package 42 in BEAST 2 v2.6.7 43 hosted on the CIPRES Science Gateway 39 . ND4 and Cytb were concatenated into a single locus to reflect the linked state of mitochondrial genes. We used the General Time Reversible + Gamma Site Model (GTR + Γ) of evolution on the data to accommodate possible rate variation among lineages with an uncorrelated log-normal relaxed clock model.
We used the Fossilized Birth-Death (FBD) model 44 for the tree prior, which provides an alternative to the incoherent approach of specifying calibration constraints and the associated probability distributions assigned to interior nodes. The FBD model directly incorporates all fossil data in a Bayesian framework 44 and requires priors for only four parameters, namely: speciation rate (λ), extinction rate (μ), fossil recovery rate (ψ), and proportion of sampled extant species (ρ). Prior values for λ and μ were adopted from previously published speciation and extinction rates for Micrurus 45 .
Based on the above four parameters the diversification rate was assigned an exponential prior with a mean value of 0.14. In order to reduce the occurrence of extreme values (i.e., 0 or 1) for the sampling proportion rate, a beta prior distribution with shape parameters α = 2.0 and β = 2.0 was assigned. Given the inherent challenge of assigning a prior for the turnover (or relative extinction) rate, we kept the default uniform prior for this parameter. The parameter ρ was set to a value of 0.58 since we had 43 of the currently considered 85 extant New World coralsnake species represented in our phylogeny. Both the hyperparameters for the clock models were assigned an exponential prior distribution, with the mean (ucldMean) set to a value of 10.0 and the standard deviation (ucldStdev) set to 1.0.
We rooted the tree prior to an origin time of 23.0 Mya based on the age of the oldest known elapid fossils 46,47 and assigned a lognormal prior distribution. Using the Mean in Real Space option the mean was specified to have a value of 8.2 and standard deviation of 1.0. The offset was set to 14.80 (i.e., the age of the oldest fossil used for calibration). Stem Micrurus were constrained with a date of 14.8 (16.0-13.6) Mya using Micrurus sp. indet, a specimen from the Barstovian of Nebraska 48,49 , which is considered the oldest New World elapid fossil. A second calibration point, a fossil dated to 830 K (1.8-0.3) Mya, was used to constrain the split between M. fulvius and M. tener as it is considered the MRCA between the two Micrurus species 50,51 . The latter age estimate has also been independently corroborated by a recent molecular study investigating the demographic history of the M. fulvius species complex 17 . Both fossil calibrations were included as blank sequences (represented by "?") in the nexus alignment file and assigned uniform distribution priors (following 44 ).
We ran two independent analyses for 1 × 10 8 generations, sampling every 1,000 generations. We used the program Tracer v2.6.7 52 to confirm the stationarity of the MCMC, with acceptable effective sample sizes of the posterior (ESS > 200) for each estimated prior. For each independent analysis, we also executed a run that samples only from priors while ignoring sequence data, which allowed us to compare marginal prior distributions on each parameter and assess the informativeness of the priors. After discarding 25% of samples as burn-in, the trees and parameter estimates from the independent runs were combined using LogCombiner v2.6.7 for a total of 10,000 final trees. Because the phylogenetic placement of the fossils was not meaningful, given that it was www.nature.com/scientificreports/ only informed by the prior (i.e., no associated sequence data), we removed them from the posterior trees using the FulltoExtantTreeConverter application in BEAST 2. Next, we used the program TreeAnnotator v2.6.7 to summarize parameter values of the samples from the posterior on a maximum clade credibility tree with mean node heights. The final tree was visualized and modified in FigTree v1.4.4 53 . The tree figures were generated in Inkscape Project (2020) (retrieved from https:// inksc ape. org).

Results
Phylogenetic relationships. This work is aligned with previous phylogenetic hypotheses of two primary sister clades of Micrurus: a triadal and a monadal clade 15,16,40,54 , with slight differences likely due to taxonomic coverage. Our phylogenetic tree recovers the triadal clade, constituted largely by South American species, and a monadal clade that includes species from South and Central America (Fig. 1).
Our results indicate that Micrurus nigrocinctus is polyphyletic due to the strongly supported inclusion of the clades constituted by M. latifasciatus nuchalis and M. ruatanus + M. n. divaricatus. All specimens of M. n. nigrocinctus from Panama (lineage 1) form a highly supported monophyletic group sister to M. mosquitensis from Costa Rica. This Isthmian Central American coralsnake group is sister to all other M. nigrocinctus ssp. represented in our work, including M. l. nuchalis. Within this group of mostly Nuclear Central American snakes, M. l. nuchalis is sister to the rest, including the well-supported clade formed by M. ruatanus + M. n. divaricatus (lineage 2) and its sister clade of all remaining M. n. nigrocinctus + M. n. zunilensis (lineage 3). The M. n. nigrocinctus + M. n. zunilensis clade is strongly supported by specimens from Southern Mexico to Costa Rica (Fig. 1).
Genetic divergence. As could be expected for intraspecific genetic differentiation, the recovered genetic divergence within each of the M. n. nigrocinctus clades was low. Divergence within the larger clade of M. n. zunilensis + M. n. nigrocintus (lineage 3) was 1.1%, within Panama M. n. nigrocintus (lineage 1) 2.2%, within M. mosquitensis 1.3%, and within the clade containing M. ruatanus 1.1%. Both M. ruatanus samples from Islas de Bahia (Roatán, Honduras) and the M. n. divaricatus skin sample (unknown locality) recovered the same haplotype. In contrast the genetic differentiation of these and M. n. divaricatus from mainland Honduras (Olancho) was 1.8% (lineage 2) (Supplementary Table S4).

Discussion
Lineage 1, Micrurus nigrocinctus from Panama. The well-supported Panamanian clade from our phylogenetic analysis strongly suggests that this group of snakes represents a separate evolutionary unit when compared to all other samples in the M. nigrocinctus species complex (Fig. 3). The oldest name available for the Panamanian group is M. nigrocinctus (Girard, 1854) 19 , a taxon with type locality in Taboga Island, a locality we have not sampled. Still, it is only 9.4 km from mainland Panama and on the Pacific discharge of the Panama Canal, in the Bay of Panama. Our molecular work and sampling efforts suggest that M. nigrocinctus is restricted to Panama, but we have a large geographical sampling gap between southwestern Costa Rica and Central Panama which might limit our conclusions. Including nuclear markers and comprehensive morphological analyses (beyond the current scope of this work) could improve our conclusion and corroborate the groups proposed by our analysis.
Micrurus ruatanus is considered a distinctive species of coralsnake, possessing a bicolored head, body and tail of alternating black and red bands (somewhat orange when small), the nuchal band not reaching parietals, and 33-45 black rings on body 64 . Micrurus ruatanus is said to be different from bicolored and tricolored mainland Honduras M. nigrocinctus (including M. n. divaricatus) with the latter having fewer black body bands (30 or fewer, according to 64 ) and females of the latter having fewer ventral scales (193-204 versus 206-224; males overlapping 183-194 versus 178-191 in M. ruatanus) 64 . Nonetheless, these differences seem minimal to separate them as different species, suggesting local differences likely conditioned through insular adaptation.
New World coralsnakes show considerable levels of introgression, as shown through ddRAD sequence data when compared to mtDNA phylogenies 11  www.nature.com/scientificreports/ taxonomy of these species raises interesting questions regarding the species´ morphology and venom evolution. The venom of Micrurus ruatanus is one of the most toxic among Central and South American snakes 65 and has been compared to the Costa Rican M. nigrocinctus proteome 65 , but further work is needed to assess its systematics and taxonomy. The M. ruatanus sequenced for this study (specimen MJ1507) is the same individual for which the venom proteome was characterized 65 . Our molecular clock calibration suggests a possible evolutionary time of less than 1.4 Mya of venom evolution from mainland M. n. divaricatus. This time frame for venom evolution could be key to examining venom toxicity in the genus and deepening our temporal understanding of venom proteome change.
The island of Roatán has several endemic reptiles, some of which also occur on nearby islets. The endemics include Anolis roatanensis, Ctenosaura oedirhina, Marisora roatanae, Norops roatanensis, Oxybelis wilsoni, Phyllodactylus palmeus, Sphaerodactylus leonardovaldesi, Sphaerodactylus rosaurae and even a mammal (the Roatán Agouti, Dasyprocta ruatanica) 66 . Stepping-stone colonization likely emerged from the Honduran islands of Cayos Cochinos and/or the island of Utila, situated on the continental shelf surrounded by shallow waters of about 30-55 m. The island chain would have been connected to the mainland during the Pleistocene, between 13,000 and 18,000 years ago, at the end of the Wisconsin glacial period 31 , indicating recent isolation 67 . Therefore, M. ruatanus might have used the same colonization routes to reach Roatán in the Pliocene-Pleistocene, when the transmarine crossing would have been as little as 20 km. Nevertheless, Roatán Island was likely never connected to the mainland. The genetic divergence between the mainland and island form suggests genetic divergence dating to no earlier than the mid-Pleistocene. The resulting over-inflation of species can be avoided by validating the molecularly delimited species through independent lines of evidence (i.e., multiple data types). Therefore, statistical-morphological and genetic-multilocus studies must delimit this species geographically and morphologically.
Micrurus mosquitensis. Micrurus nigrocinctus mosquitensis was described from Limón, Costa Rica, by Schmidt in 1933 59 and treated as a subspecies for almost 70 years, until it was given full species status in 2004 by Solorzano 26 . Our phylogenetic analyses, including Pacific and Atlantic populations, suggest this taxon as distinct and the Pacific and Atlantic versant populations as belonging to the same taxon (Fig. 3). Using molecular data we confirm the arrangement 26 , contrary to McCranie 64 .
Divergence estimates. The timings suggest that the New World coralsnakes emerged at the onset of the Langhian period. Micrurus dates to circa 16 Mya, when both monadal and triadal clades split. The triadal clade dates circa 15 Mya while the younger monadal clade likely dates to circa 10 Mya (Fig. 2). Overall, these fossilbased time estimates to infer timing for the whole genus remain congruent to other studies based on fossil estimates 16,45,[68][69][70] in comparison to when they are inferred with other dating methods 10,71 . Although our time estimates for the triadal clade are similar to those of Lee et al. 45 and Zaher et al. 68,69 , they remain discordant to those of Hurtado Gomez et al. 70 . In addition, in their study 70 they recover a more ancestral branching of the monadal rather than the triadal clade, in contradiction with other published works (e.g., 69 ). Their non-inclusion  www.nature.com/scientificreports/ of the triadal snake M. laticollaris in their study likely influenced their tree topology. Similarly, the monadal clade dating remains older in our analyses than in previous studies, though similar to two previous works 10,45 .
The divergence between species within the M. nigrocinctus group falls within the Late Pliocene and Early Pleistocene, suggesting a time of climatic and environmental changes favorable to speciation events in the region. The low mtDNA genetic divergence between Roatán M. ruatanus and northern Honduras M. n. divaricatus suggests recent connections between localities, most likely through a recent Pleistocene dispersal event. Until this study, only M. diutius had been sequenced from an island (Trinidad) and compared to their counterpart in Guyana, recovering no genetic divergence for the ND4 locus 10 . This study further increases our understanding of colonization processes of New World coralsnakes on islands and in this instance, points to a recent connection due to a combination of low glacial sea levels during the Pleistocene and island hopping across deep (~ 500 m), but narrow (~ 20-40 km) ocean channels.
It must be noted that a certain degree of discrepancy is expected between results generated by different studies, given the differences in the sampling of genes, taxa, calibration types (fossil vs secondary), and notably the choice of priors for key parameters in BEAST. Additionally, node ages estimated under the FBD model are typically incongruent with those resulting from traditional node dating analyses as consistently shown by several recent studies 45,72,73 . Biomedical implications. Most of the published works on the composition and toxicity of the venom in Micrurus nigrocinctus come from specimens from the Central Plateau of Costa Rica (e.g., 3,20,74,75 ). In contrast, studies on the toxicity and composition of the venom of M. nigrocinctus in Panama are lacking. However, some preliminary chromatographic profiles suggest a high similarity with venoms of M. n. nigrocinctus from Costa Rica (data not shown) and some preliminary data from the research center and information on medicines and toxins (Centro de Investigación e Información en Medicamentos y Tóxicos, personal communication) of Panama confirm patient recovery from M. nigrocintus bites with Costa Rican antivenom. Fortunately, the most extensive geographic range within the M. nigrocinctus group is occupied by our lineage 3, which also occupies Central Costa Rica, suggesting that all of the antivenoms produced could have applicability for bites throughout Central America. Micrurus venoms can be classified as either PLA2-rich or 3FTx-rich 74,76 and cross-species antivenom efficacy is highly dependent on the compositional predominance of either. Micrurus nigrocinctus from lineage 3, the most important coralsnakes from the medical point of view in Central America, has a venom rich in PLA2 77 and the anticoral antivenom manufactured in Costa Rica uses this venom as an immunogen. This makes this antivenom mainly effective against snake venoms rich in PLA2, although it is less effective against those rich in 3FTx 1,78 . For example, both the venom of M. mosquitensis, a Costa Rican species sympatric with M. n. nigrocinctus lineage 3, predominantly rich in PLA2 8 , as well as that of M. ruatanus, predominantly rich in 3FTx, are neutralized by the antivenom produced in Costa Rica 65 .
Current limitations and further taxonomical implications. The findings of this study are based on mitochondrial derived data and future research on these clades would greatly benefit from the addition of informative nuclear data. Unfortunately, nuclear markers for this study proved problematic to sequence and those that amplified were too conservative to resolve with confidence phylogenetic inferences. Therefore, the nuclear data were not included in the final analyses to avoid the improper statement that the phylogeny was derived from both genomes. Fortunately, a new ML tree RAD-seq data (Renan J. Bosque personal communication) derived from some of the exact same Micrurus individuals used for the mitochondrial DNA tree of this study is fully congruent with the recovered phylogenetic positions of lineages 1, 2 and 3, and therefore suggests that the mitochondrial tree is a species tree. Future phylogenetic congruence (may it be from morphological, -RAD-seq or nuclear sequence data) will suggest the need for taxonomical changes of these lineages.

Conclusion.
Our results strongly highlight divergent mtDNA clades within the M. nigrocinctus species complex. The evidence supported by morphological and/or nuclear gene-based evidence will likely elevate them to species level. Micrurus nigrocinctus from Panama is monophyletic and highly divergent from all other lineages. Using genetics we confirm the already elevated species M. mosquitensis 26 in Costa Rica. Micrurus ruatanus seems restricted to Honduras and closely related to M. n. divaricatus from the mainland, with very low genetic differentiation. A group of currently recognized M. n. nigrocinctus and M. n zulinensis are recovered as monophyletic, with low genetic divergence despite their broad distribution from southern Mexico to Costa Rica, including Guatemala, Honduras, and Nicaragua. Furthermore, the phylogenetic placement of the two samples of M. latifasciatus nuchalis within the M. nigrocinctus species complex confirms the paraphyly of the group and therefore the need for a taxonomic rearrangement. Our findings have pivotal implications for future antivenom research, but we refrain from performing taxonomic changes until more data can be assessed for this complex of taxa.

Data availability
The datasets generated and/or analyzed during the current study are available in the GenBank repository (https:// www. ncbi. nlm. nih. gov/ genba nk/).